Index: /tags/photdbc-1-0/Ohana/src/photdbc/Makefile =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/Makefile (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/Makefile (revision 2617) @@ -0,0 +1,74 @@ +include ../../Configure +HOME = $(ROOT)/src/photdbc +CONFIG = $(ROOT)/config +PROGRAM = photdbc + +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) -lohana -lFITS -lm +CFLAGS = -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +SRCLIST = \ +$(SRC)/photdbc.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/name_region.$(ARCH).o \ +$(SRC)/find_images.$(ARCH).o \ +$(SRC)/find_regions.$(ARCH).o \ +$(SRC)/get_regions.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/misc.$(ARCH).o \ +$(SRC)/time.$(ARCH).o \ +$(SRC)/join_stars.$(ARCH).o \ +$(SRC)/unique_measures.$(ARCH).o \ +$(SRC)/flag_measures.$(ARCH).o \ +$(SRC)/get_mags.$(ARCH).o \ +$(SRC)/check_directory.$(ARCH).o \ +$(SRC)/subcatalog.$(ARCH).o \ +$(SRC)/wcatalog.$(ARCH).o + +OBJ = $(SRCLIST) + +default: $(PROGRAM) + +$(SRCLIST): $(INC)/photdbc.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/photdbc-1-0/Ohana/src/photdbc/bin/.cvsignore =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/bin/.cvsignore (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/bin/.cvsignore (revision 2617) @@ -0,0 +1,1 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix Index: /tags/photdbc-1-0/Ohana/src/photdbc/doc/notes.txt =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/doc/notes.txt (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/doc/notes.txt (revision 2617) @@ -0,0 +1,28 @@ + +- select a region + +- load average / measure / secfilt / missing + +- construct subset based on selection criteria: + - Minst + - dMinst + - X,Y + - R,D + +- merge neighbor stars (join within X.X arcsec) + +- mark bad measurements: + + - outside area + - Minst < cutoff + - Mcal > 0.015 (lookup from table?) + + - subset for given filter (match by string?): + - 2x 3sigma + - softened chisq > 3.0 + +- mark objects with fewer than N measurements + + +----- + Index: /tags/photdbc-1-0/Ohana/src/photdbc/include/photdbc.h =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/include/photdbc.h (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/include/photdbc.h (revision 2617) @@ -0,0 +1,103 @@ +# include +# include +# include + +typedef struct { + int measure; // pointer to measure entry + unsigned int averef; // old measures come from multiple averages + float R, D; // actual coordinates of this measure +} Mpointer; + +typedef struct { + double median; + double mean; + double sigma; + double error; + double chisq; + double min; + double max; + double total; + int Nmeas; +} StatType; + +// need to get RADIUS from Config + +/* global variables */ +int SHOW_PARAMS; +int VERBOSE; +double JOIN_RADIUS; +double UNIQ_RADIUS; +double DMCAL_MIN;; +char ImageCat[256]; +char GSCFILE[256]; +char CATDIR[256]; +char PhotCodeFile[256]; + +double RMIN; +double RMAX; +double DMIN; +double DMAX; + +double XMIN; +double XMAX; +double YMIN; +double YMAX; +double MMIN; +double MMAX; + +double DMSYS; +double DMGAIN; +double CHISQ_MAX; +double SIGMA_MAX; +int NMEAS_MIN; + +PhotCodeData photcodes; + +# define FLAG_AREA 0X0001 +# define FLAG_MINST 0X0002 +# define FLAG_DMCAL 0x0004 +# define FLAG_SIGCLIP 0X0008 +# define FLAG_CHISQ 0x0010 +# define FLAG_TOOFEW 0x0020 +# define FLAG_DUPMEAS 0x0040 + +Image *find_images (GSCRegion *region, int Nregion, int *Nimages, int **LineNum); +GSCRegion *find_regions (Image *image, int Nimage, int *Nregions, GSCRegion *fullregion); +Image *get_images (int *N); +GSCRegion *get_regions (double minRa, double maxRa, double minDec, double maxDec, int *Nregions); +Catalog *load_catalogs (GSCRegion *region, int Nregion); +GSCRegion *load_images (char *seed, int *nregion); +GSCRegion *name_region (char *name, int *Nregions); +char *sec_to_date (unsigned long second); +void ConfigInit (int *argc, char **argv); +int Fseek (FILE *f, long offset, int whence); +int args (int argc, char **argv); +int chk_time (char *line); +int corner_check (double *x1, double *y1, double *x2, double *y2); +int dms_to_ddd (double *Value, char *string); +int edge_check (double *x1, double *y1, double *x2, double *y2); +void flag_measures (Catalog *catalog); +int gcatalog (Catalog *catalog); +void get_mags (Catalog *catalog); +void getfullregion (Image *image, int Nimage, GSCRegion *fullregion); +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string); +void initialize (int argc, char **argv); +void initstats (char *mode); +void join_stars (Catalog *catalog); +int liststats (double *value, double *dvalue, int N, StatType *stats); +unsigned long jd_to_sec (double jd); +unsigned long date_to_sec (char *date); +double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); +double sec_to_jd (unsigned long second); +void set_ZP (double ZERO); +void sortA (double *X, int N); +void sortB (double *X, double *Y, int N); +void sortC (double *X, double *Y, double *F1, double *F2, int N); +void sortD (double *X, double *Y, double *Z, int N); +void sort_lists (float *X, float *Y, int *index, int N); +void sort_time (unsigned int *X, int *Y, int N); +int str_to_dtime (char *line, double *second); +int str_to_radec (double *ra, double *dec, char *str1, char *str2); +int str_to_time (char *line, unsigned int *second); +void unique_measures (Catalog *catalog); +void wcatalog (Catalog *catalog); Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/ConfigInit.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/ConfigInit.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/ConfigInit.c (revision 2617) @@ -0,0 +1,71 @@ +# include "photdbc.h" + +int success = TRUE; + +void WarnConfig (char *config, char *key, char * mode, int N, void *var) { + if (!ScanConfig (config, key, mode, N, var)) { + fprintf (stderr, "missing config variable %s\n", key); + success = FALSE; + } +} + +void ConfigInit (int *argc, char **argv) { + + double ZERO_POINT, tmp; + 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); + + WarnConfig (config, "JOIN_RADIUS", "%lf", 0, &JOIN_RADIUS); + WarnConfig (config, "UNIQ_RADIUS", "%lf", 0, &UNIQ_RADIUS); + + WarnConfig (config, "XMIN", "%lf", 0, &XMIN); + WarnConfig (config, "XMAX", "%lf", 0, &XMAX); + WarnConfig (config, "YMIN", "%lf", 0, &YMIN); + WarnConfig (config, "YMAX", "%lf", 0, &YMAX); + + WarnConfig (config, "RMIN", "%lf", 0, &RMIN); + WarnConfig (config, "RMAX", "%lf", 0, &RMAX); + WarnConfig (config, "DMIN", "%lf", 0, &DMIN); + WarnConfig (config, "DMAX", "%lf", 0, &DMAX); + + WarnConfig (config, "DMCAL_MIN", "%lf", 0, &tmp); DMCAL_MIN = 1000*tmp; + WarnConfig (config, "MMIN", "%lf", 0, &tmp); MMIN = 1000*tmp; + WarnConfig (config, "MMAX", "%lf", 0, &tmp); MMAX = 1000*tmp; + WarnConfig (config, "DMSYS", "%lf", 0, &tmp); DMSYS = SQ(1000*tmp); + WarnConfig (config, "DMGAIN", "%lf", 0, &DMGAIN); + WarnConfig (config, "SIGMA_MAX", "%lf", 0, &SIGMA_MAX); SIGMA_MAX *= 1000; + + WarnConfig (config, "CHISQ_MAX", "%lf", 0, &CHISQ_MAX); + WarnConfig (config, "NMEAS_MIN", "%d", 0, &NMEAS_MIN); + + WarnConfig (config, "IMAGE_CATALOG", "%s", 0, ImageCat); + WarnConfig (config, "GSCFILE", "%s", 0, GSCFILE); + WarnConfig (config, "CATDIR", "%s", 0, CATDIR); + WarnConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + WarnConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + + if (!success) { + fprintf (stderr, "missing config parameter\n"); + exit (1); + } + + SetZeroPoint (ZERO_POINT); + + if (!LoadPhotcodes (PhotCodeFile)) { + fprintf (stderr, "error loading photcodes\n"); + exit (1); + } + + free (config); + free (file); + +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/args.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/args.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/args.c (revision 2617) @@ -0,0 +1,25 @@ +# include "photdbc.h" + +int args (int argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + SHOW_PARAMS = FALSE; + if ((N = get_argument (argc, argv, "-params"))) { + remove_argument (N, &argc, argv); + SHOW_PARAMS = TRUE; + } + + if (argc != 2) { + fprintf (stderr, "USAGE: photdbc (region)\n"); + exit (2); + } + +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/check_directory.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/check_directory.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/check_directory.c (revision 2617) @@ -0,0 +1,64 @@ +# include "photdbc.h" + +check_directory (char *basefile) { + + FILE *f; + char *c, dir[256], filename[256]; + struct stat filestat; + uid_t fuid, uid; + gid_t fgid, gid; + int status, cmode; + + 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, "directory %s does not exist, creating...\n", dir); + cmode = S_IRWXU | S_IRWXG | S_IRWXO; + status = mkdir (dir, cmode); + if (status == -1) { + fprintf (stderr, "ERROR: can't create %s\n", dir); + exit (1); + } + } + status = stat (dir, &filestat); + 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 (1); + } +} + +/* uses cp only */ +make_backup (char *filename) { + + int status, cmode; + struct stat filestat; + char line [256]; + + status = stat (filename, &filestat); + if (status == 0) { /* file exists, make backup copy */ + sprintf (line, "cp %s %s~\0", filename, filename); + status = system (line); + if (status) { + fprintf (stderr, "ERROR: unable to create %s~, exiting\n", filename); + exit (1); + } + cmode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + sprintf (line, "%s~\0", filename); + chmod (line, cmode); + } + +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/find_images.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/find_images.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/find_images.c (revision 2617) @@ -0,0 +1,238 @@ +# include "photdbc.h" + +/* this function returns a list of all images which overlap the given + set of region files. All images in the image catalog are tested + once, so there is no check that an image already has been included. + LineNum stores the locations in the Image database of the list of + images */ + +Image *find_images (GSCRegion *region, int Nregion, int *Nimages, int **LineNum) { + + Header header; + Image *timage, *image; + int i, j, k, m, found, nimage, Nimage, NIMAGE, NTIMAGE, Nloop, Nlast; + int n, Nim, status, InRange, np, Np; + FILE *f; + double Xc[5], Yc[5], Xi[5], Yi[5], r, d, dx, dy; + char line[256]; + int *line_number; + Coords tcoords; + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + /* we make positional comparisons in the projection of catalog */ + tcoords.crval1 = 0.5*(region[0].RA[0] + region[0].RA[1]); + tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]); + tcoords.crpix1 = tcoords.crpix2 = 0.0; + 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"); + + /* 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); + + /* load all images */ + ALLOCATE (timage, Image, Nimage); + status = Fread (timage, sizeof(Image), Nimage, f, "image"); + if (status != Nimage) { + fprintf (stderr, "ERROR: couldn't read images from image catalog: %s\n", ImageCat); + exit (0); + } + + nimage = 0; + NIMAGE = 100; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, int, NIMAGE); + + for (i = 0; i < Nimage; i++) { + +# if (0) + /* select images by photcode */ + np = timage[i].source; + Np = photcodes.hashcode[np]; + if (Np < 0) continue; + if (photcodes.code[Np].equiv != PHOTCODE) continue; + + /* select images by time */ + if (TimeSelect) { + if (timage[i].tzero < TSTART) continue; + if (timage[i].tzero > TSTOP) continue; + } +# endif + + /* 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; + found = FALSE; + /* transform to tcoords */ + for (j = 0; j < 5; j++) { + XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords); + InRange = RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords); + if (!InRange) goto imskip; + } + /* compare with each region file */ + for (m = 0; (m < Nregion) && !found; m++) { + /* define catalog corners */ + Xc[0] = region[m].RA[0]; Yc[0] = region[m].DEC[0]; + Xc[1] = region[m].RA[1]; Yc[1] = region[m].DEC[0]; + Xc[2] = region[m].RA[1]; Yc[2] = region[m].DEC[1]; + Xc[3] = region[m].RA[0]; Yc[3] = region[m].DEC[1]; + Xc[4] = region[m].RA[0]; Yc[4] = region[m].DEC[0]; + for (j = 0; j < 5; j++) { + r = Xc[j]; d = Yc[j]; + RD_to_XY (&Xc[j], &Yc[j], r, d, &tcoords); + } + 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; + + /* check if image corner inside catalog */ + for (j = 0; (j < 4) && !found; j++) { + found |= corner_check (&Xi[j], &Yi[j], &Xc[0], &Yc[0]); + } + /* check if catalog corner inside image */ + for (j = 0; (j < 4) && !found; j++) { + found |= corner_check (&Xc[j], &Yc[j], &Xi[0], &Yi[0]); + } + /* check if edges cross */ + for (j = 0; (j < 4) && !found; j++) { + for (k = 0; (k < 4) && !found; k++) { + found |= edge_check (&Xi[j], &Yi[j], &Xc[k], &Yc[k]); + } + } + if (found) { + image[nimage] = timage[i]; + if (image[nimage].code == ID_IMAGE_NOCAL) { + image[nimage].code &= ~ID_IMAGE_NOCAL; + } + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, int, NIMAGE); + } + } + } + imskip: + continue; + } + + if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); + + /* we are going to use image[nimage] to represent images not in the database: + these are measurements from external sources, like USNO */ + + assignMcal (&image[nimage], (double *) NULL, -1); + image[nimage].code = ID_IMAGE_NEW; + + REALLOCATE (image, Image, MAX (nimage + 1, 1)); + REALLOCATE (line_number, int, MAX (nimage + 1, 1)); + free (timage); + *Nimages = nimage; + *LineNum = line_number; + fclose (f); + return (image); +} + +/* check if line between points 0 and 1 of x1 + crosses line between points 0 and 1 of x2 */ +int edge_check (double *x1, double *y1, double *x2, double *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); + +} + +/* check if point x1,y1 is in box formed by x2[0-4] */ +int corner_check (double *x1, double *y1, double *x2, double *y2) { + + int i; + double theta; + + theta = 0; + + for (i = 0; i < 4; i++) { + theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); + } + if (fabs(theta) > 6) { + return (TRUE); + } else { + return (FALSE); + } +} + +/* returns the opening angle between the three points (2 is in middle) + in range -pi to pi */ + +double opening_angle (double x1, double y1, double x2, double y2, double x3, double 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); + +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/find_regions.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/find_regions.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/find_regions.c (revision 2617) @@ -0,0 +1,102 @@ +# include "photdbc.h" + +/* given a list of images, find all region files which overlap these images. + if a region file does not exist, we should recognize this fact somehow */ + +GSCRegion *find_regions (Image *image, int Nimage, int *Nregions, GSCRegion *fullregion) { + + int i, nregions; + double x, y; + GSCRegion *regions; + double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec; + + if (VERBOSE) fprintf (stderr, "finding regions: "); + + /* find RA and DEC range of images */ + MinRA = MinDEC = +361; + MaxRA = MaxDEC = -361; + for (i = 0; i < Nimage; i++) { + 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[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); + MinRA = MIN (ra, MinRA); + MaxRA = MAX (ra, MaxRA); + MinDEC = MIN (dec, MinDEC); + MaxDEC = MAX (dec, MaxDEC); + } + } + } + + /* if images overlap 0,360 boundary, then use -180 to +180 range */ + if (MaxRA - MinRA > 180) { + MinRA = +361; + MaxRA = -361; + for (i = 0; i < Nimage; i++) { + 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[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); + if (ra > 180) ra -= 360; + MinRA = MIN (ra, MinRA); + MaxRA = MAX (ra, MaxRA); + } + } + } + } + + fullregion[0].RA[0] = MinRA; + fullregion[0].RA[1] = MaxRA; + fullregion[0].DEC[0] = MinDEC; + fullregion[0].DEC[1] = MaxDEC; + + /* find all region files containing specified ra, dec range */ + regions = get_regions (MinRA, MaxRA, MinDEC, MaxDEC, &nregions); + + if (VERBOSE) fprintf (stderr, "%d found\n", nregions); + + *Nregions = nregions; + return (regions); +} + + +void getfullregion (Image *image, int Nimage, GSCRegion *fullregion) { + + int i; + double x, y; + double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec; + + /* find RA and DEC range of images */ + MinRA = MinDEC = +361; + MaxRA = MaxDEC = -361; + for (i = 0; i < Nimage; i++) { + 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[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); + MinRA = MIN (ra, MinRA); + MaxRA = MAX (ra, MaxRA); + MinDEC = MIN (dec, MinDEC); + MaxDEC = MAX (dec, MaxDEC); + } + } + } + + /* if images overlap 0,360 boundary, then use -180 to +180 range */ + if (MaxRA - MinRA > 180) { + MinRA = +361; + MaxRA = -361; + for (i = 0; i < Nimage; i++) { + 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[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); + if (ra > 180) ra -= 360; + MinRA = MIN (ra, MinRA); + MaxRA = MAX (ra, MaxRA); + } + } + } + } + + fullregion[0].RA[0] = MinRA; + fullregion[0].RA[1] = MaxRA; + fullregion[0].DEC[0] = MinDEC; + fullregion[0].DEC[1] = MaxDEC; +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/flag_measures.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/flag_measures.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/flag_measures.c (revision 2617) @@ -0,0 +1,203 @@ +# include "photdbc.h" + +void flag_measures (Catalog *catalog) { + + int i, j, k, m, n, N, found; + int *imagelink, Nimage, Nmeasure, Nlist, PhotNsec; + unsigned int Time; + int Ncode, Photcode, Minst, Mag, BAD_MEASURE; + double ra, dec, x, y, Nsigma; + double *list, *dlist; + + Image *image; + StatType stats; + + if (VERBOSE) fprintf (stderr, "flagging bad measurements\n"); + + /* load images */ + image = get_images (&Nimage); + + PhotNsec = GetPhotcodeNsecfilt (); + /* need to replace dMcal with values from external table? */ + + /* clear existing measure flags */ + for (i = 0; i < catalog[0].Nmeasure; i++) catalog[0].measure[i].flags = 0; + + /* generate image index */ + ALLOCATE (imagelink, int, catalog[0].Nmeasure); + for (i = 0; i < catalog[0].Nmeasure; i++) { + n = catalog[0].measure[i].averef; + Time = catalog[0].measure[i].t; + Photcode = catalog[0].measure[i].source; + found = FALSE; + for (k = 0; !found && (k < Nimage); k++) { + if ((Time == image[k].tzero) && (Photcode == image[k].source)) { + imagelink[i] = k; + found = TRUE; + } + } + if (!found) { + /* some error? */ + } + } + + /* flag by area : find source image, get X,Y coords, set flag */ + for (i = 0; i < catalog[0].Nmeasure; i++) { + n = catalog[0].measure[i].averef; + ra = catalog[0].average[n].R - catalog[0].measure[i].dR / 360000.0; + dec = catalog[0].average[n].D - catalog[0].measure[i].dD / 360000.0; + k = imagelink[i]; + RD_to_XY (&x, &y, ra, dec, &image[k].coords); + if (x < XMIN) goto flagarea; + if (x > XMAX) goto flagarea; + if (y < YMIN) goto flagarea; + if (y > YMAX) goto flagarea; + continue; + flagarea: + catalog[0].measure[i].flags |= FLAG_AREA; + } + + /* flag by Minst */ + for (i = 0; i < catalog[0].Nmeasure; i++) { + Minst = iPhotInst (&catalog[0].measure[i]); + if (Minst < MMIN) goto flagminst; + if (Minst > MMAX) goto flagminst; + continue; + flagminst: + catalog[0].measure[i].flags |= FLAG_MINST; + } + + /* flag by image (dMcal? time from table?) */ + for (i = 0; i < catalog[0].Nmeasure; i++) { + k = imagelink[i]; + if (image[k].dMcal > DMCAL_MIN) { + catalog[0].measure[i].flags |= FLAG_DMCAL; + } + } + + /* allocate a list for temp storage of mag values */ + Nlist = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + Nlist = MAX (Nlist, catalog[0].average[i].Nm); + } + ALLOCATE (list, double, Nlist); + ALLOCATE (dlist, double, Nlist); + initstats ("MEAN"); + + /* the following sections apply only to PRI/SEC photcodes */ + for (n = 0; n < photcodes.Ncode; n++) { + if (photcodes.code[n].type == PHOT_DEP) continue; + if (photcodes.code[n].type == PHOT_REF) continue; + Photcode = photcodes.code[n].code; + + BAD_MEASURE = FLAG_AREA | FLAG_MINST | FLAG_DMCAL; + /* 2x 3sigma clipping */ + for (i = 0; i < catalog[0].Naverage; i++) { + + /* extract good measures, calculate stats */ + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + if (catalog[0].measure[m].flags & BAD_MEASURE) continue; + list[N] = iPhotRel (&catalog[0].measure[m], &catalog[0].average[i], &catalog[0].secfilt[i*PhotNsec]); + dlist[N] = sqrt (SQ(DMGAIN*catalog[0].measure[m].dM) + DMSYS); + N++; + } + if (N < 3) continue; + liststats (list, dlist, N, &stats); + Nsigma = 3.0; + if (N < 10) Nsigma = 0.9*N / sqrt (N - 1); + + /* re-extract good measures, calculate stats */ + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + if (catalog[0].measure[m].flags & BAD_MEASURE) continue; + Mag = iPhotRel (&catalog[0].measure[m], &catalog[0].average[i], &catalog[0].secfilt[i*PhotNsec]); + if (fabs(Mag - stats.mean) > Nsigma*stats.sigma) continue; + list[N] = Mag; + dlist[N] = sqrt (SQ(DMGAIN*catalog[0].measure[m].dM) + DMSYS); + N++; + } + if (N > 2) liststats (list, dlist, N, &stats); + + /* flag bad measures */ + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + Mag = iPhotRel (&catalog[0].measure[m], &catalog[0].average[i], &catalog[0].secfilt[i*PhotNsec]); + if (fabs(Mag - stats.mean) > Nsigma*stats.sigma) catalog[0].measure[m].flags |= FLAG_SIGCLIP; + } + } + + BAD_MEASURE = FLAG_AREA | FLAG_MINST | FLAG_DMCAL | FLAG_SIGCLIP; + /* flag variables */ + for (i = 0; i < catalog[0].Naverage; i++) { + /* extract good measures, calculate chisq */ + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + if (catalog[0].measure[m].flags & BAD_MEASURE) continue; + list[N] = iPhotRel (&catalog[0].measure[m], &catalog[0].average[i], &catalog[0].secfilt[i*PhotNsec]); + dlist[N] = sqrt (SQ(DMGAIN*catalog[0].measure[m].dM) + DMSYS); + N++; + } + liststats (list, dlist, N, &stats); + if (stats.chisq < CHISQ_MAX) continue; + + /* flag bad measures */ + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + catalog[0].measure[m].flags |= FLAG_CHISQ; + } + } + + /* flag too few valid measures */ + BAD_MEASURE = FLAG_AREA | FLAG_MINST | FLAG_DMCAL | FLAG_SIGCLIP | FLAG_CHISQ; + for (i = 0; i < catalog[0].Naverage; i++) { + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + if (catalog[0].measure[m].flags & BAD_MEASURE) continue; + N++; + } + if (N == 0) continue; // only flag objects with valid measurements, but too few + if (N >= NMEAS_MIN) continue; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + catalog[0].measure[m].flags |= FLAG_TOOFEW; + } + } + } +} + + +/* notes + + flag by area: uses an inefficient search for the image. + images should be sorted by time and we should use bisection + + number of loops over Nmeasure: 9 + +*/ + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/gcatalog.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/gcatalog.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/gcatalog.c (revision 2617) @@ -0,0 +1,22 @@ +# include "photdbc.h" + +int gcatalog (Catalog *catalog) { + + char mode, status; + + mode = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (!load_catalog (catalog, mode, VERBOSE)) { + fprintf (stderr, "ERROR: failure loading catalog %s\n", catalog[0].filename); + exit (1); + } + + if (VERBOSE) { + fprintf (stderr, "read %d stars from catalog file %s (%d measurements)\n", + catalog[0].Naverage, catalog[0].filename, catalog[0].Nmeasure); + } + return (TRUE); +} + + + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/get_mags.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/get_mags.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/get_mags.c (revision 2617) @@ -0,0 +1,65 @@ +# include "photdbc.h" + +void get_mags (Catalog *catalog) { + + int i, j, m, n, N; + int Nlist, Nsec, Nsecfilt, Ncode; + short int *Mval; + int Photcode, BAD_MEASURE; + double *list, *dlist; + StatType stats; + + if (VERBOSE) fprintf (stderr, "re-calculating photometry\n"); + + BAD_MEASURE = FLAG_AREA | FLAG_MINST | FLAG_SIGCLIP | FLAG_CHISQ | FLAG_TOOFEW; + Nsecfilt = catalog[0].Nsecfilt; + + /* allocate a list for temp storage of mag values */ + Nlist = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + Nlist = MAX (Nlist, catalog[0].average[i].Nm); + catalog[0].average[i].M = NO_MAG; + catalog[0].average[i].dM = NO_MAG; + catalog[0].average[i].Xm = NO_MAG; + for (j = 0; j < Nsecfilt; j++) { + catalog[0].secfilt[i*Nsecfilt+j].M = NO_MAG; + catalog[0].secfilt[i*Nsecfilt+j].dM = NO_MAG; + catalog[0].secfilt[i*Nsecfilt+j].Xm = NO_MAG; + } + } + ALLOCATE (list, double, Nlist); + ALLOCATE (dlist, double, Nlist); + initstats ("WT_MEAN"); + + /* check on photcode */ + for (n = 0; n < photcodes.Ncode; n++) { + if (photcodes.code[n].type == PHOT_DEP) continue; + if (photcodes.code[n].type == PHOT_REF) continue; + Photcode = photcodes.code[n].code; + Nsec = photcodes.hashNsec[Photcode]; + + /* find average magnitudes */ + for (i = 0; i < catalog[0].Naverage; i++) { + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + if ((Ncode = photcodes.hashcode[catalog[0].measure[m].source]) < 0) continue; + if (photcodes.code[Ncode].equiv != Photcode) continue; + if (photcodes.code[Ncode].type != PHOT_DEP) continue; + if (catalog[0].measure[m].flags & BAD_MEASURE) continue; + list[N] = PhotRel (&catalog[0].measure[m], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); + dlist[N] = sqrt (SQ(DMGAIN*catalog[0].measure[m].dM) + DMSYS); + N++; + } + liststats (list, dlist, N, &stats); + if (N < 1) continue; + + Mval = (Nsec == -1) ? &catalog[0].average[i].M : &catalog[0].secfilt[i*Nsecfilt+Nsec].M; + *Mval = stats.mean; + Mval = (Nsec == -1) ? &catalog[0].average[i].dM : &catalog[0].secfilt[i*Nsecfilt+Nsec].dM; + *Mval = stats.sigma; + Mval = (Nsec == -1) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[i*Nsecfilt+Nsec].Xm; + *Mval = 100.0*log10(stats.chisq); + } + } +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/get_regions.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/get_regions.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/get_regions.c (revision 2617) @@ -0,0 +1,208 @@ +# include "photdbc.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}; + +double Areas[] = {2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0, + 2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.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}; + +/* returns a list of region files within the desired RA, DEC region */ +GSCRegion *get_regions (double minRa, double maxRa, double minDec, double maxDec, int *Nregions) { + + char buffer[28800], temp[50]; + GSCRegion *regions; + FILE *f; + double RA0, RA1, DEC0, DEC1; + int i, j, NBigDec; + int NLINES, done, NREGIONS, nregions; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find regions file %s\n", GSCFILE); + exit (0); + } + + NREGIONS = 50; + ALLOCATE (regions, GSCRegion, NREGIONS); + nregions = 0; + + if ((minDec <= -90) || (maxDec >= 90)) { + minRa = 0; + maxRa = 360; + } + + /* use the pole regions, if near pole */ + if (maxDec > 86.25) { + sprintf (regions[nregions].filename, "n8230/pole.cpt\0"); + regions[nregions].RA[0] = 0; + regions[nregions].RA[1] = 360; + regions[nregions].DEC[0] = 86.25; + regions[nregions].DEC[1] = 90.0; + nregions ++; + if (nregions == NREGIONS) { + NREGIONS += 50; + REALLOCATE (regions, GSCRegion, NREGIONS); + } + } + + if (minDec > 86.25) { + *Nregions = nregions; + return (regions); + } + + if ((minDec < 0) && (maxDec > 0)) { + /* Search Both Sides */ + NBigDec = 0; + } else { + /* find large DEC region (directory) */ + NBigDec = -1; + for (i = 0; i < 12; i++) { + if ((minDec >= BigDecBounds[i]) && (minDec < BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + if (NBigDec < 0) { + for (i = 13; i < 24; i++) { + if ((maxDec < BigDecBounds[i]) && (maxDec >= BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + } + } + if (NBigDec < 0) { + fprintf (stderr, "ERROR: Dec out of range: %f\n", minDec); + *Nregions = 0; + return ((GSCRegion *) NULL); + } + + /* count lines before section */ + NLINES = 0; + for (i = 0; i < NBigDec; i++) { + NLINES += NDecLines[i]; + } + fseek (f, 5*2880 + 48*NLINES, SEEK_SET); + + /* should be in this section. if not, there is a problem counting... */ + /* careful with the 0,360.0 boundary **/ + done = FALSE; + for (j = 0; !done && (NBigDec + j < 25); j++) { + fread (buffer, 48*NDecLines[NBigDec + j], 1, f); + for (i = 0; (i < NDecLines[NBigDec + j]); 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; + if ((DEC1 > 0) && (minDec < DEC1) && (maxDec > DEC0) && (minRa < RA1) && (maxRa > RA0)) { + temp[5] = 0; + sprintf (regions[nregions].filename, "%s/%s.cpt\0", Dec2Sections[NBigDec + j], &temp[1]); + regions[nregions].RA[0] = RA0; + regions[nregions].RA[1] = RA1; + regions[nregions].DEC[0] = DEC0; + regions[nregions].DEC[1] = DEC1; + nregions ++; + if (nregions == NREGIONS) { + NREGIONS += 50; + REALLOCATE (regions, GSCRegion, NREGIONS); + } + } + if ((DEC1 < 0) && (minDec < DEC0) && (maxDec > DEC1) && (minRa < RA1) && (maxRa > RA0)) { + temp[5] = 0; + sprintf (regions[nregions].filename, "%s/%s.cpt\0", Dec2Sections[NBigDec + j], &temp[1]); + regions[nregions].RA[0] = RA0; + regions[nregions].RA[1] = RA1; + regions[nregions].DEC[0] = DEC0; + regions[nregions].DEC[1] = DEC1; + nregions ++; + if (nregions == NREGIONS) { + NREGIONS += 50; + REALLOCATE (regions, GSCRegion, NREGIONS); + } + } + if (((DEC1 > 0) && (maxDec <= DEC1)) || ((DEC1 < 0) && (minDec >= DEC1))) { + done = TRUE; + } + } + if (done && (minDec < 0) && (maxDec > 0) && (BigDecBounds[NBigDec + j + 1] > 0)) { + /* skip remaining north sections, try south sections */ + /* count lines before section */ + NLINES = 0; + for (i = 0; i < 13; i++) { + NLINES += NDecLines[i]; + } + fseek (f, 5*2880 + 48*NLINES, SEEK_SET); + done = FALSE; + j = 12; + } + } + + REALLOCATE (regions, GSCRegion, MAX (1,nregions)); + *Nregions = nregions; + + fclose (f); + return (regions); + +} + +/**********/ +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) { + + 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); +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/initialize.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/initialize.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/initialize.c (revision 2617) @@ -0,0 +1,38 @@ +# include "photdbc.h" + +void initialize (int argc, char **argv) { + + int Np; + double tmp; + + /* are these set correctly? */ + ConfigInit (&argc, argv); + args (argc, argv); + + if (SHOW_PARAMS) { + fprintf (stderr, "current parameter settings:\n"); + fprintf (stderr, "VERBOSE: %d\n", VERBOSE); + fprintf (stderr, "JOIN_RADIUS: %lf\n", JOIN_RADIUS); + fprintf (stderr, "UNIQ_RADIUS: %lf\n", UNIQ_RADIUS); + + fprintf (stderr, "XMIN: %lf\n", XMIN); + fprintf (stderr, "XMAX: %lf\n", XMAX); + fprintf (stderr, "YMIN: %lf\n", YMIN); + fprintf (stderr, "YMAX: %lf\n", YMAX); + fprintf (stderr, "MMIN: %lf\n", MMIN); + fprintf (stderr, "MMAX: %lf\n", MMAX); + fprintf (stderr, "DMCAL_MIN: %lf\n", DMCAL_MIN); + fprintf (stderr, "DMSYS: %lf\n", DMSYS); + + fprintf (stderr, "CHISQ_MAX: %lf\n", CHISQ_MAX); + fprintf (stderr, "NMEAS_MIN: %lf\n", NMEAS_MIN); + + fprintf (stderr, "IMAGE_CATALOG %s\n", ImageCat); + fprintf (stderr, "GSCFILE %s\n", GSCFILE); + fprintf (stderr, "CATDIR %s\n", CATDIR); + fprintf (stderr, "PHOTCODE_FILE %s\n", PhotCodeFile); + + exit (0); + } +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/join_stars.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/join_stars.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/join_stars.c (revision 2617) @@ -0,0 +1,184 @@ +# include "photdbc.h" + +void join_stars (Catalog *catalog) { + + int i, j, k, m, M, Ni, Nj, first_j, Ns1, Ns2; + int Naves, Nmeas, Ncurr; + int Naverage, Nmeasure, *found, *index; + float *X, *Y, dX, dY, dR, RADIUS2; + double Sr, Sd, Rmid, Dmid; + + Average *naverage, *average; + Measure *nmeasure, *measure; + SecFilt *nsecfilt, *secfilt; + Mpointer *mpointer; + Coords tcoords; + + if (VERBOSE) fprintf (stderr, "joining overlapping stars\n"); + + average = catalog[0].average; + measure = catalog[0].measure; + secfilt = catalog[0].secfilt; + + Naverage = catalog[0].Naverage; + Nmeasure = catalog[0].Nmeasure; + RADIUS2 = JOIN_RADIUS*JOIN_RADIUS; + + /* reference for coords is this region */ + Rmid = Dmid = 0; + for (i = 0; i < Naverage; i++) { + Rmid += average[i].R; + Dmid += average[i].D; + } + Rmid /= Naverage; + Dmid /= Naverage; + + /* coordinate system for projection */ + tcoords.crval1 = Rmid; + tcoords.crval2 = Dmid; + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 0; + strcpy (tcoords.ctype, "RA---TAN"); + + /* project & sort coordinates in local linear frame */ + ALLOCATE (X, float, Naverage); + ALLOCATE (Y, float, Naverage); + ALLOCATE (index, int, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + fRD_to_XY (&X[i], &Y[i], average[i].R, average[i].D, &tcoords); + } + sort_lists (X, Y, index, Naverage); + + /* flags to mark if a star has been handled */ + ALLOCATE (found, int, Naverage); + for (i = 0; i < Naverage; i++) found[i] = FALSE; + + ALLOCATE (naverage, Average, Naverage); + ALLOCATE (mpointer, Mpointer, Nmeasure); + for (i = 0; i < Nmeasure; i++) mpointer[i].averef = -1; + + Naves = 0; // counter for new averages + Nmeas = 0; // counter for new measures + for (i = j = 0; (i < Naverage) && (j < Naverage); ) { + + Ni = index[i]; + Nj = index[j]; + + /* a new star, add it to naverage[] */ + if (!found[i]) { + naverage[Naves] = average[Ni]; + naverage[Naves].offset = Nmeas; + for (k = 0; k < average[Ni].Nm; k++) { + m = average[Ni].offset + k; + mpointer[Nmeas].measure = m; + mpointer[Nmeas].averef = Naves; + mpointer[Nmeas].R = average[Ni].R - measure[m].dR / 360000.0; + mpointer[Nmeas].D = average[Ni].D - measure[m].dD / 360000.0; + Nmeas ++; + } + Ncurr = Naves; + found[i] = TRUE; + Naves ++; + } + + if (found[j]) { j++; continue; } // don't duplicate + if (j == i) { j++; continue; } // don't auto-correlate + + dX = X[i] - X[j]; + if (dX <= -JOIN_RADIUS) { /* X[j] is too large */ + while (found[i] && (i < Naverage)) i++; + continue; + } + if (dX >= +JOIN_RADIUS) { /* X[i] is too large */ + j++; + continue; + } + + first_j = j; + for (; (dX > -2*JOIN_RADIUS) && (j < Naverage); j++) { + Nj = index[j]; + if (found[j]) continue; + dX = X[i] - X[j]; + dY = Y[i] - Y[j]; + dR = dX*dX + dY*dY; + if (dR < RADIUS2) { /* matched star, join to first */ + + /* define pointers for new measures */ + for (k = 0; k < average[Nj].Nm; k++) { + m = average[Nj].offset + k; + mpointer[Nmeas].measure = m; + mpointer[Nmeas].averef = Ncurr; + mpointer[Nmeas].R = average[Nj].R - measure[m].dR / 360000.0; + mpointer[Nmeas].D = average[Nj].D - measure[m].dD / 360000.0; + Nmeas ++; + } + naverage[Ncurr].Nm += average[Nj].Nm; + found[j] = TRUE; + + /* recalculate naverage[Ncurr].RA,DEC */ + /* this must be done here to keep the average position consistent + for the next star found */ + for (Sr = Sd = k = 0; k < naverage[Ncurr].Nm; k++) { + m = naverage[Ncurr].offset + k; + Sr += mpointer[m].R; + Sd += mpointer[m].D; + } + Sr = Sr / naverage[Ncurr].Nm; + Sd = Sd / naverage[Ncurr].Nm; + naverage[Ncurr].R = Sr; + naverage[Ncurr].D = Sd; + + /* update original measurement offsets */ + for (k = 0; k < naverage[Ncurr].Nm; k++) { + m = naverage[Ncurr].offset + k; + M = mpointer[m].measure; + measure[M].dR = 360000.0*(Sr - mpointer[m].R); + measure[M].dD = 360000.0*(Sd - mpointer[m].D); + } + + /* update current reference star position */ + fRD_to_XY (&X[i], &Y[i], Sr, Sd, &tcoords); + } + } + while (found[i] && (i < Naverage)) i++; + j = first_j; + } + if (Nmeas != Nmeasure) { + fprintf (stderr, "failure to match measures\n"); + } + + /* create a new Measure table in the appropriate sequence */ + ALLOCATE (nmeasure, Measure, Nmeasure); + for (i = 0; i < Nmeasure; i++) { + nmeasure[i] = measure[mpointer[i].measure]; + nmeasure[i].averef = mpointer[i].averef; + } + + /* create an empty SecFilt table */ + ALLOCATE (nsecfilt, SecFilt, Naves*catalog[0].Nsecfilt); + bzero (nsecfilt, Naves*catalog[0].Nsecfilt*sizeof(SecFilt)); + + free (catalog[0].average); + free (catalog[0].measure); + free (catalog[0].secfilt); + + catalog[0].average = naverage; + catalog[0].measure = nmeasure; + catalog[0].secfilt = nsecfilt; + + catalog[0].Naverage = Naves; +} + +/* notes: + merge averages of stars within RADIUS to single stars + updates the average RA and DEC fields + creates a new Average table 'naverage' + updates the Measure table to the new sequence + generates a new SecFilt table with empty values + does NOT update magnitudes +*/ + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/liststats.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/liststats.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/liststats.c (revision 2617) @@ -0,0 +1,114 @@ +# include "photdbc.h" + +enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, + M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN}; + +static int statmode; + +void initstats (char *mode) { + + statmode = -1; + if (!strcmp (mode, "MEAN")) statmode = M_MEAN; + if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN; + if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN; + if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN; + if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN; + if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN; + if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN; + + if (statmode == -1) { + fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode); + exit (1); + } +} + +int liststats (double *value, double *dvalue, int N, StatType *stats) { + + int i, ks, ke, Nm; + double Mo, dMo, M, dM, X2, dS, *chi; + + stats[0].Nmeas = N; + stats[0].mean = 0; + if (N < 2) return (FALSE); + + sortB (value, dvalue, N); + stats[0].median = value[(int)(0.5*N)]; + stats[0].min = value[0]; + stats[0].max = value[N-1]; + + switch (statmode) { + case M_MEDIAN: + ks = 0; + ke = N; + Mo = stats[0].median; + Nm = N; + goto chisq; + break; + case M_MEAN: + case M_WT_MEAN: + ks = 0; + ke = N; + break; + case M_INNER_MEAN: + case M_INNER_WTMEAN: + case M_CHI_INNER_MEAN: + case M_CHI_INNER_WTMEAN: + ks = 0.25*N + 0.50; + ke = 0.75*N + 0.25; + if (N <= 3) { + ks = 0; + ke = N; + } + break; + } + + if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) { + ALLOCATE (chi, double, N); + for (i = 0; i < N; i++) { + chi[i] = (value[i] - stats[0].median) / dvalue[i]; + } + sortD (chi, value, dvalue, N); + free (chi); + } + + /* calculating the per-star offset based on the weighted average */ + M = dM = Nm = 0; + if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) { + for (i = ks; i < ke; i++) { + M += value[i] / SQ (dvalue[i]); + dM += 1.0 / SQ (dvalue[i]); + Nm ++; + } + Mo = M / dM; + dMo = sqrt (1.0 / dM); + } else { + for (i = ks; i < ke; i++) { + M += value[i]; + dM += SQ (dvalue[i]); + Nm ++; + } + Mo = M / (double) Nm; + dMo = sqrt (dM / (double) Nm); + } + + chisq: + /* find sigma and chisq */ + X2 = dS = 0; + for (i = ks; i < ke; i++) { + M = SQ (value[i] - Mo); + dM = SQ (dvalue[i]); + X2 += M / dM; + dS += M; + } + X2 = X2 / Nm; + dS = sqrt (dS / Nm); + + stats[0].mean = Mo; + stats[0].Nmeas = Nm; + stats[0].chisq = X2; + stats[0].sigma = dS; + stats[0].error = dMo; + + return (TRUE); +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/load_catalogs.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/load_catalogs.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/load_catalogs.c (revision 2617) @@ -0,0 +1,33 @@ +# include "photdbc.h" + +Catalog *load_catalogs (GSCRegion *region, int Nregion) { + + int i, Nstar; + Catalog *catalog; + + if (VERBOSE) fprintf (stderr, "loading catalog data\n"); + + ALLOCATE (catalog, Catalog, Nregion); + + /* load data from each region file */ + for (i = 0; i < Nregion; i++) { + catalog[i].filename = region[i].filename; /* don't free region before catalog! */ + switch (lock_catalog (&catalog[i], LCK_SOFT)) { + case 0: + fprintf (stderr, "ERROR: can't lock file %s\n", catalog[i].filename); + exit (1); + case 1: + gcatalog (&catalog[i]); /* load from disk */ + break; + case 2: + if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog[i].filename); + catalog[i].Naverage = 0; + catalog[i].Nmeasure = 0; + break; + default: + fprintf (stderr, "weird lock_catalog exit state\n"); + exit (1); + } + } + return (catalog); +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/load_images.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/load_images.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/load_images.c (revision 2617) @@ -0,0 +1,64 @@ +# include "photdbc.h" + +static int Nimage; +static int *Line; +static Image *image; + +GSCRegion *load_images (char *seed, int *nregion) { + + GSCRegion *region; + int Nregion, Ntimage; + Image *timage; + + region = name_region (seed, &Nregion); + + /* find images which supplied all measurements */ + image = find_images (region, Nregion, &Nimage, &Line); + *nregion = Nregion; + + return (region); +} + +Image *get_images (int *N) { + + *N = Nimage; + return (image); +} + +/* +void wimages () { + + int i, status, mode; + FILE *f; + Header header; + + setMcalFinal (image, Nimage); + + if (!fits_read_header (ImageCat, &header)) { + fprintf (stderr, "ERROR: can't find %s\n", ImageCat); + exit (0); + } + + f = fopen (ImageCat, "r+"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't open image catalog file to append: %s\n", ImageCat); + exit (0); + } + mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + chmod (ImageCat, mode); + + if (VERBOSE) fprintf (stderr, "writing out %d images\n", Nimage); + for (i = 0; i < Nimage; i++) { + if (image[i].code == ID_IMAGE_NOCAL) image[i].code &= ~ID_IMAGE_NOCAL; + + fseek (f, header.size + Line[i]*sizeof(Image), SEEK_SET); + status = Fwrite (&image[i], sizeof(Image), 1, f, "image"); + if (status != 1) { + fprintf (stderr, "ERROR: failure writing image data to catalog %s\n", ImageCat); + exit (0); + } + } + + fclose (f); +} +*/ Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/misc.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/misc.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/misc.c (revision 2617) @@ -0,0 +1,273 @@ + +void sortA (double *X, int N) { + + int l,j,ir,i; + double tX; + + if (N < 2) return; + + 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; + } +} + +void sortB (double *X, double *Y, int N) { + + int l,j,ir,i; + double tX, tY; + + if (N < 2) return; + + 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; + } +} + +void sortC (double *X, double *Y, double *F1, double *F2, int N) { + + int l,j,ir,i; + double tX, tY; + double t1, t2; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + t1 = F1[l]; + t2 = F2[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + t1 = F1[ir]; + F1[ir] = F1[0]; + t2 = F2[ir]; + F2[ir] = F2[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + F1[0] = t1; + F2[0] = t2; + 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]; + F1[i] = F1[j]; + F2[i] = F2[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + F1[i] = t1; + F2[i] = t2; + } +} + +void sortD (double *X, double *Y, double *Z, int N) { + + int l,j,ir,i; + double tX, tY, tZ; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tZ = Z[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tZ = Z[ir]; + Z[ir] = Z[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + Z[0] = tZ; + 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]; + Z[i] = Z[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + Z[i] = tZ; + } +} + +void sort_lists (float *X, float *Y, int *S, int N) { + + int l,j,ir,i; + double tX, tY; + int tS; + + if (N < 2) return; + + 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; + } +} + +void sort_time (unsigned int *X, int *Y, int N) { + + int l,j,ir,i; + unsigned int tX; + int tY; + + if (N < 2) return; + + 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; + } +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/name_region.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/name_region.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/name_region.c (revision 2617) @@ -0,0 +1,24 @@ +# include "photdbc.h" + +GSCRegion *name_region (char *name, int *Nregions) { + + GSCRegion *region; + Header header; + + ALLOCATE (region, GSCRegion, 1); + + sprintf (region[0].filename, "%s/%s", CATDIR, name); + if (!fits_read_header (region[0].filename, &header)) { + if (VERBOSE) fprintf (stderr, "catalog file not yet available: %s\n", region[0].filename); + exit (0); + } + fits_scan (&header, "RA0", "%lf", 1, ®ion[0].RA[0]); + fits_scan (&header, "RA1", "%lf", 1, ®ion[0].RA[1]); + fits_scan (&header, "DEC0", "%lf", 1, ®ion[0].DEC[0]); + fits_scan (&header, "DEC1", "%lf", 1, ®ion[0].DEC[1]); + fits_free_header (&header); + + *Nregions = 1; + return (region); +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/photdbc.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/photdbc.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/photdbc.c (revision 2617) @@ -0,0 +1,52 @@ +# include "photdbc.h" + +int main (int argc, char **argv) { + + int Nregion; + char filename[1024]; + Catalog *catalog; + GSCRegion *region; + + /* get configuration info, args, lockfile */ + initialize (argc, argv); + + /* load region and images overlapping region */ + region = load_images (argv[1], &Nregion); + + /* load catalog data corresponding to region file */ + catalog = load_catalogs (region, Nregion); + + /* limit number of measures based on selections */ + catalog = subcatalog (catalog); + + /* limit number of measures based on selections */ + join_stars (catalog); + + unique_measures (catalog); + + flag_measures (catalog); + + get_mags (catalog); + + unlock_catalog (catalog); + + sprintf (filename, "newdir/%s", argv[1]); + catalog[0].filename = filename; + + check_directory (catalog[0].filename); + if (!lock_catalog (catalog, LCK_XCLD)) { + fprintf (stderr, "can't lock new catalog file %s\n", catalog[0].filename); + exit (1); + } + wcatalog (catalog); +} + + +/** this program is somewhat rough: minimal error checking + - + only one catalog is loaded, corresponding to region file + + + + + */ Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/subcatalog.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/subcatalog.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/subcatalog.c (revision 2617) @@ -0,0 +1,107 @@ +# include "photdbc.h" + +Catalog *subcatalog (Catalog *catalog) { + + int i, j, offset, ecode; + int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; + int mag; + Catalog *subcatalog; + + ALLOCATE (subcatalog, Catalog, 1); + + Nsecfilt = GetPhotcodeNsecfilt (); + + subcatalog[0] = catalog[0]; + + /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ + NAVERAGE = 50; + NMEASURE = 1000; + ALLOCATE (subcatalog[0].average, Average, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + ALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + Nmeasure = Naverage = 0; + + for (i = 0; i < catalog[0].Naverage; i++) { + /* exclude stars with too few measurements */ + if (catalog[0].average[i].Nm < NMEAS_MIN) continue; + + /* exclude stars not in range */ + if (catalog[0].average[i].R < RMIN) continue; + if (catalog[0].average[i].R > RMAX) continue; + if (catalog[0].average[i].D < DMIN) continue; + if (catalog[0].average[i].D > DMAX) continue; + + /* assign average and secfilt values */ + subcatalog[0].average[Naverage] = catalog[0].average[i]; + subcatalog[0].average[Naverage].offset = Nmeasure; + for (j = 0; j < Nsecfilt; j++) { + subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j]; + } + + Nm = 0; + for (j = 0; j < catalog[0].average[i].Nm; j++) { + + offset = catalog[0].average[i].offset + j; + + /* exclude measurements by photcode + if (DropPhotcode) { + ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source); + if (ecode == photcode[0].code) continue; + }*/ + + /* select measurements by quality + if (DropNonStellar && catalog[0].measure[offset].dophot != 1) continue; + */ + + /* exclude measurements by measurement error */ + if (catalog[0].measure[offset].dM > SIGMA_MAX) continue; + + /* select measurements by mag limit */ + mag = iPhotInst (&catalog[0].measure[offset]); + if (mag < MMIN) continue; + if (mag > MMAX) continue; + + subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; + subcatalog[0].measure[Nmeasure].averef = Naverage; + Nmeasure ++; + Nm ++; + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + } + } + if (Nm < NMEAS_MIN) { /* enough measurements meet criteria? */ + Nmeasure -= Nm; + continue; + } + subcatalog[0].average[Naverage].Nn = 0; + subcatalog[0].average[Naverage].Nm = Nm; + Naverage ++; + if (Naverage == NAVERAGE) { + NAVERAGE += 50; + REALLOCATE (subcatalog[0].average, Average, NAVERAGE); + REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + } + } + REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1)); + subcatalog[0].Naverage = Naverage; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + + ALLOCATE (subcatalog[0].missing, Missing, 1); + subcatalog[0].Nmissing = 0; + + if (catalog[0].Naverage > 0) free (catalog[0].average); + if (catalog[0].Nmeasure > 0) free (catalog[0].measure); + if (catalog[0].Naverage > 0) free (catalog[0].secfilt); + + if (VERBOSE) { + fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, + subcatalog[0].Naverage, subcatalog[0].Nmeasure); + } + + return (subcatalog); + +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/time.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/time.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/time.c (revision 2617) @@ -0,0 +1,313 @@ +# include + +/***** convert [-]00:00:00 to 0.0000 ****/ +int dms_to_ddd (double *Value, char *string) { + + int valid, neg, status; + double tmp, value; + char *p1, *p2, *px; + + valid = FALSE; + neg = FALSE; + stripwhite (string); + p1 = string; + px = string + strlen(string); + + if (string[0] == '-') { + valid = TRUE; + neg = TRUE; + p1 = &string[1]; + } + if (string[0] == '+') { + valid = TRUE; + neg = FALSE; + p1 = &string[1]; + } + if (isdigit(string[0])) { + valid = TRUE; + p1 = &string[0]; + } + if (!valid) { return (FALSE); } + + status = 1; + tmp = strtod (p1, &p2); + if (p2 == p1) return (FALSE); /* entry not a number: +fred */ + value = tmp; + if (p2 == px) goto escape; /* entry only number: +1.0 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:fred */ + status = 2; + value += tmp / 60.0; + if (p2 == px) goto escape; /* entry only number: +1:1 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */ + value += tmp / 3600.0; + + escape: + if (neg) { + value *= -1; + } + *Value = value; + + return (status); +} + +/**********/ +int str_to_radec (double *ra, double *dec, char *str1, char *str2) { + + double Ra, Dec; + + *ra = *dec = 0; + switch (dms_to_ddd (&Ra, str1)) { + case 0: + fprintf (stderr, "syntax error in RA\n"); + return (FALSE); + case 1: + break; + case 2: + Ra = Ra * 15; + break; + } + switch (dms_to_ddd (&Dec, str2)) { + case 0: + fprintf (stderr, "syntax error in DEC\n"); + return (FALSE); + case 1: + case 2: + break; + } + *ra = Ra; + *dec = Dec; + return (TRUE); +} + +# define FORMAT_DAYS 1 +# define FORMAT_HOURS 2 +# define FORMAT_MINUTES 3 +# define FORMAT_SECONDS 4 +# define FORMAT_JD 5 +# define FORMAT_DATE 6 + +/**********/ +int chk_time (char *line) { + + char *p1, *p2; + double tmp; + int mode; + + p1 = line; + tmp = strtod (p1, &p2); + if (p2 == p1 + strlen (p1) - 1) { + if (*p2 == 'd') { + mode = FORMAT_DAYS; + } + if (*p2 == 'h') { + mode = FORMAT_HOURS; + } + if (*p2 == 'm') { + mode = FORMAT_MINUTES; + } + if (*p2 == 's') { + mode = FORMAT_SECONDS; + } + if (*p2 == 'j') { + mode = FORMAT_JD; + } + } else { + mode = FORMAT_DATE; + } + return (mode); +} + +/**********/ +int str_to_time (char *line, unsigned int *second) { + + struct timeval now; + double jd; + + if (!strcasecmp (line, "NOW")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = now.tv_sec; + return (TRUE); + } + + if (!strcasecmp (line, "TODAY")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = 86400 * ((int)(now.tv_sec / 86400)); + return (TRUE); + } + + switch (chk_time (line)) { + case 0: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + case FORMAT_JD: + jd = strtod (line, 0); + *second = jd_to_sec (jd); + return (TRUE); + case FORMAT_DATE: + *second = date_to_sec (line); + return (TRUE); + } + return (FALSE); +} + +/**********/ +int str_to_dtime (char *line, double *second) { + + switch (chk_time (line)) { + case 0: + case FORMAT_JD: + case FORMAT_DATE: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + } + return (FALSE); +} + +/**********/ +double sec_to_jd (unsigned long second) { + + double jd; + + jd = second/86400.0 + 2440587.5; + return (jd); +} + +/**********/ +unsigned long jd_to_sec (double jd) { + + unsigned long int second; + + second = (jd - 2440587.5)*86400; + return (second); +} + +/**********/ +char *sec_to_date (unsigned long second) { + + struct tm *gmt; + char *line; + + ALLOCATE (line, char, 64); + gmt = gmtime (&second); + sprintf (line, "%4d/%02d/%02d,%02d:%02d:%02d", 1900+gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec); + return (line); + +} + +/***** date in format yyyy/mm/dd,hh:mm:ss *****/ +unsigned long date_to_sec (char *date) { + + unsigned long second; + double tmp, jd; + struct tm now; + char *p1, *p2, *px; + + p1 = date; + px = date + strlen(date); + bzero (&now, sizeof(now)); + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_year = tmp; + if (now.tm_year > 1000) now.tm_year -= 1900; + if (now.tm_year < 50) now.tm_year += 100; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mon = tmp - 1; /* mon runs from 0 - 11 */ + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mday = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + p1 = p2 + 1; + now.tm_hour = tmp; + if (p2 == px) goto escape; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_min = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_sec = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + escape: + jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4) + + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12) + - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5; + + second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec; + + return (second); +} + + +/* times may be in forms as: + 20040200450s (N seconds since 1970.0) + 2440900.232j (julian date) + 99/02/23,03:22:18 (date string) + (separators may be anything except space, +, -) + 99:02:15:12:23:30 + 99:02:15:12h23m30s + */ + + +/* fseek with timeout - 0.5 sec */ +int Fseek (FILE *f, long offset, int whence) { + + int status, k; + + status = fseek (f, offset, whence); + if (status == -1) { + int k; + /* fprintf (stderr, "problem seeking position: %d\n", errno); */ + for (k = 0; (k < 10) && ((status = fseek (f, 0, SEEK_SET)) == -1); k++) usleep (50000); + if (status == -1) { + /* fprintf (stderr, "ERROR: serious problem seeking position: %d\n", errno); */ + /* clear lock at this point? */ + return (0); + } + } + return (1); +} Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/unique_measures.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/unique_measures.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/unique_measures.c (revision 2617) @@ -0,0 +1,101 @@ +# include "photdbc.h" + +void unique_measures (Catalog *catalog) { + + int i, j, k, m, N; + int No, Nm, Nvalid, Nlist; + int *ilist; + unsigned int *tlist; + Measure *mlist, *nmeasure; + double dR, dD, R; + Average *p; + + if (VERBOSE) fprintf (stderr, "removing duplicate measurements\n"); + + /* allocate a list for temp storage of measure values */ + Nlist = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + Nlist = MAX (Nlist, catalog[0].average[i].Nm); + } + ALLOCATE (ilist, int, Nlist); + ALLOCATE (tlist, unsigned int, Nlist); + ALLOCATE (mlist, Measure, Nlist); + + p = catalog[0].average; + for (i = 0; i < catalog[0].Naverage; i++) { + + if (catalog[0].average != p) { + fprintf (stderr, "error\n"); + } + + /* extract list of measurements */ + N = 0; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + ilist[N] = m; + tlist[N] = catalog[0].measure[m].t; + N++; + } + + /* sort the images by time */ + sort_time (tlist, ilist, N); + + /* look for duplicates */ + for (j = k = 0; (j < N) && (k < N); j++) { + if (j == k) { + k++; + continue; + } + if (tlist[j] != tlist[k]) { + j++; + k = j + 1; + continue; + } + if (catalog[0].measure[ilist[j]].flags & FLAG_DUPMEAS) { + j++; + k = j + 1; + continue; + } + + dR = (catalog[0].measure[ilist[j]].dR - catalog[0].measure[ilist[k]].dR); + dD = (catalog[0].measure[ilist[j]].dD - catalog[0].measure[ilist[k]].dD); + R = SQ(dR) + SQ(dD); + + /* dR, dD in 1/100 arcsec, RADIUS is 0.1 arcsec */ + if (R < UNIQ_RADIUS) { + /* matched pair */ + catalog[0].measure[ilist[k]].flags |= FLAG_DUPMEAS; + } + k++; + } + } + + Nvalid = 0; + for (i = 0; i < catalog[0].Nmeasure; i++) { + if (!(catalog[0].measure[i].flags & FLAG_DUPMEAS)) Nvalid ++; + } + + No = 0; + ALLOCATE (nmeasure, Measure, MAX (Nvalid, 1)); + for (i = 0; i < catalog[0].Naverage; i++) { + /* copy valid measures to new table, update offset, Nm */ + Nm = 0; + m = catalog[0].average[i].offset; + catalog[0].average[i].offset = No; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++, No++) { + if (catalog[0].measure[m].flags & FLAG_DUPMEAS) continue; + nmeasure[No] = catalog[0].measure[m]; + if (catalog[0].measure[m].averef != i) { + fprintf (stderr, "average reference mismatch!!\n"); + } + Nm ++; + } + catalog[0].average[i].Nm = Nm; + } + + free (catalog[0].measure); + catalog[0].measure = nmeasure; + catalog[0].Nmeasure = No; + +} + Index: /tags/photdbc-1-0/Ohana/src/photdbc/src/wcatalog.c =================================================================== --- /tags/photdbc-1-0/Ohana/src/photdbc/src/wcatalog.c (revision 2617) +++ /tags/photdbc-1-0/Ohana/src/photdbc/src/wcatalog.c (revision 2617) @@ -0,0 +1,25 @@ +# include "photdbc.h" + +void wcatalog (Catalog *catalog) { + + fits_modify (&catalog[0].header, "NRPHOT", "%t", 1, TRUE); + + save_catalog (catalog, VERBOSE); + + if (catalog[0].Naverage) { + free (catalog[0].average); + catalog[0].Naverage = 0; + } + if (catalog[0].Nmeasure) { + free (catalog[0].measure); + catalog[0].Nmeasure = 0; + } + if (catalog[0].Nmissing) { + free (catalog[0].missing); + catalog[0].Nmissing = 0; + } + if (catalog[0].Nsecfilt) { + free (catalog[0].secfilt); + catalog[0].Nsecfilt = 0; + } +}